function [pathwayStatistics,listOfPoolMetabolites,zeroSolution, ...
    listOfPathways,listOfAllWeights, modelToReturn]=...
        get_pathway_statistics(pathway, model, shouldModelBeReturned, isForProduct, originalFlux)

     defines = evalin('base', 'defines'); 

   if ~isForProduct
        objectiveReactionIndex = 1;
        skipObjectiveReactionFluxAnalysis = false;
        optimizationObjective = 'max';
        if str2num(defines.isSubstrateUsed)
           objectiveReactionIndex = str2num(defines.substrateObjectiveReaction);
        elseif str2num(defines.isProductUsed)
           objectiveReactionIndex = str2num(defines.productObjectiveReaction);
        end

        if objectiveReactionIndex == 1
           skipObjectiveReactionFluxAnalysis = true;
        end
   end
    
   if str2num(defines.isProductUsed) & str2num(defines.isSubstrateUsed)
        skipObjectiveReactionFluxAnalysis = true;
   end
    
    if strcmp(defines.isThermodinamicUsed, 'true')
        thermodinamicsValues = evalin('base', 'thermodinamicsValues'); 
    end

    elementsToExclude = '';
    if evalin('base','exist(''elementsToExclude'',''var'')')
        elementsToExclude = evalin('base','elementsToExclude');
    end
    
    currentModel = model;
    allMetabolitesInPathway = [];
    addedMetabolites = [];
    existingMetabolites = [];
    listOfPoolMetabolites = [];
    sourceTransportReaction ='';
    zeroSolution = false;
    source_metabolite = pathway(1,:);
    listOfPathways = {};
    modelToReturn = '';
    
    pathwayStatistics = struct( ...
                            'MET_COUNT',0, ...
                            'EXIST_MET_COUNT',0, ...
                            'ADDED_MET_COUNT',0, ...
                            'REAC_COUNT',0, ...
                            'EXISTING_REACT_COUNT',0, ...
                            'REACTIONS_WITH_KNOWN_THERMODINAMICS',0, ...
                            'UNBALANCED_REACT_COUNT',0, ...
                            'FBA_RESULT', 0, ... 
                            'FBA_RESULT_OBJECTIVE_REACTION', 0, ...  
                            'FBA_BIOMASS_RESULT', 0 );

    for y = 1:length(pathway(:,1))
       if regexp(pathway(y,:), 'R') == 1
           pathway(y,:)
            listOfAnalyzedMetabolites = [];
            pathwayStatistics.REAC_COUNT = pathwayStatistics.REAC_COUNT + 1;
            reactionStructure = get_reaction_information(pathway(y,:))
 
            for j = 1:length(reactionStructure.PATHWAY)
                if length(listOfPathways) > 0
                    [isReactionMemeberOfPathways, index] = ismember(reactionStructure.PATHWAY{j}, listOfPathways(:,1));
                else
                    isReactionMemeberOfPathways = false;
                end

                if isReactionMemeberOfPathways == 0
                   listOfPathways = [listOfPathways; {reactionStructure.PATHWAY{j}}, 1];
                else
                   listOfPathways(index,2) = {listOfPathways{index,2} + 1};
                end
            end

            equation = getfield(reactionStructure, 'EQUATION');              
            metabolites = regexp(equation, 'C.....', 'match');

            metabolitesCountInReaction = 0;
            metaboliteCountAddedToModel = 0;
            metabolitesCountExistingInModel = 0;

            for n = 1:size(metabolites,2)
                metaboliteStructure = struct;

               [metaboliteAbbreviation,isFoundInModel] =  ...
                   get_metabolite_abbreviation_from_model(currentModel, metabolites{n}, defines.cystoleCompartment);

               temporaryMetabolite = metabolites{n};
               
                if isempty(strmatch(metabolites{n}, pathway))
                    allMetabolitesInPathway = [allMetabolitesInPathway; metabolites{n}];
                end
                if length(metaboliteAbbreviation) > 12
                    if regexp(equation, metaboliteAbbreviation(1:12))
                        temporaryMetabolite = metaboliteAbbreviation;
                    end
                end
                equation = strrep(equation, temporaryMetabolite, metaboliteAbbreviation);
                metabolitesCountInReaction = metabolitesCountInReaction + 1;
                if isFoundInModel == 1
                    if isempty(~strmatch(metabolites{n},existingMetabolites))
                       if isempty(regexp(metaboliteAbbreviation, 'added', 'once'))
                           metabolitesCountExistingInModel = metabolitesCountExistingInModel + 1;           
                           existingMetabolites = [existingMetabolites; metabolites{n}];
                       end
                    end
                else
                    if ~strcmp(metabolites{n},elementsToExclude)
                        addedMetabolites = [addedMetabolites; metabolites{n}];
                        metaboliteCountAddedToModel = metaboliteCountAddedToModel + 1;
                    end
                end
                metaboliteStructure = get_compound_information(metabolites{n})
                listOfAnalyzedMetabolites = [listOfAnalyzedMetabolites; {metaboliteAbbreviation, metaboliteStructure}];
            end
            
            reactionStructure = setfield(reactionStructure,'EQUATION_FROM_MODEL', equation);
            reactionStructure = setfield(reactionStructure,'metabolitesCountInReaction', metabolitesCountInReaction);
            reactionStructure = setfield(reactionStructure,'metabolitesCountExistingInModel', metabolitesCountExistingInModel);
            reactionStructure = setfield(reactionStructure,'metaboliteCountAddedToModel', metaboliteCountAddedToModel);

            name = getfield(reactionStructure, 'NAME');
            reactionStructure.NAME = regexprep(name,' ','_');
            pathwayStatistics.MET_COUNT = pathwayStatistics.MET_COUNT + metabolitesCountInReaction;
            pathwayStatistics.EXIST_MET_COUNT = pathwayStatistics.EXIST_MET_COUNT + metabolitesCountExistingInModel;
            pathwayStatistics.ADDED_MET_COUNT = pathwayStatistics.ADDED_MET_COUNT + metaboliteCountAddedToModel;

            currentReactionEquation = equation;
            reactionsNumberWithThermodinamics = 0;
            
            if strcmp(defines.isThermodinamicUsed,'true')
                isThermodinamicsValueFound = is_thermodinamics_data_found(pathway(y,:));
                reactionsNumberWithThermodinamics = reactionsNumberWithThermodinamics + isThermodinamicsValueFound;
            end
            
            reactionName = reactionStructure.NAME
            [currentModel, iFExist] = addReaction(currentModel, {pathway(y,:),reactionName}, currentReactionEquation); 
                
            if iFExist > 0
                pathwayStatistics.EXISTING_REACT_COUNT = pathwayStatistics.EXISTING_REACT_COUNT + 1;
            end

            for z = 1:length(listOfAnalyzedMetabolites)
                if regexp(listOfAnalyzedMetabolites{z,1}, 'added_')
                    if regexp(listOfAnalyzedMetabolites{z,1}, source_metabolite)
                        tmpTransport = strcat('R_added_', source_metabolite,defines.cystoleCompartment, '_transport');
                        if ismember(tmpTransport, currentModel.rxns)
                            sourceTransportReaction = tmpTransport;
                        else
                            sourceTransportReaction = strcat('R_', source_metabolite, '_transport')
                        end
                        if ~ismember(sourceTransportReaction, currentModel.rxns)                     
                           exctracellularSourceMetabolite = strcat('added_',source_metabolite, defines.extracellularCompartment);                   
                           transportReactionEquation = sprintf('%s %s %s',listOfAnalyzedMetabolites{z,1}, '->', exctracellularSourceMetabolite)
                           
            
                           currentModel = addReaction(currentModel, sourceTransportReaction, transportReactionEquation); 

                           sourceMetaboliteExchangeReaction = strcat('EX_',source_metabolite,'(e)');
                           exchangeReactionEquation = sprintf('%s %s %s',exctracellularSourceMetabolite, '<=>');
                           currentModel = addReaction(currentModel, sourceMetaboliteExchangeReaction, exchangeReactionEquation); 

                           sourceMetaboliteIndex = find(strcmp(listOfAnalyzedMetabolites{z,1}, currentModel.mets),1);
                           sourceExtracellularMetaboliteIndex = find(strcmp(exctracellularSourceMetabolite, currentModel.mets),1);

                           currentMetaboliteStructure = listOfAnalyzedMetabolites(z, 2);

                           currentModel.metFormulas{sourceMetaboliteIndex} = currentMetaboliteStructure{1}.FORMULA;
                           currentModel.metFormulas{sourceExtracellularMetaboliteIndex} = currentMetaboliteStructure{1}.FORMULA;

                           keggID = regexp(currentMetaboliteStructure{1}.ENTRY, 'C.....', 'match');
                           currentModel.metKEGGID{sourceMetaboliteIndex} = keggID{1};
                           currentModel.metKEGGID{sourceExtracellularMetaboliteIndex} = keggID{1};
                        end
                    else
                        pool_metabolite_transport_rxn = strcat('R_', listOfAnalyzedMetabolites{z,1},'_transport');
                        if ~ismember(pool_metabolite_transport_rxn, currentModel.rxns)  
                            isMetaboliteFoundInPathway = false;
                            for x = 1:length(pathway(:,1))
                                if regexp( listOfAnalyzedMetabolites{z,1}, pathway(x,:))>0
                                   isMetaboliteFoundInPathway = true;
                                end
                            end
                            currentMetaboliteStructure = listOfAnalyzedMetabolites(z, 2);
                            if ~isMetaboliteFoundInPathway
                                extracellulatrPoolMetabolite = strcat(strrep(listOfAnalyzedMetabolites{z,1}, defines.cystoleCompartment, defines.extracellularCompartment));
                                transportReactionEquation = sprintf('%s %s %s',listOfAnalyzedMetabolites{z,1}, '<=>', extracellulatrPoolMetabolite);
                                currentModel = addReaction(currentModel, pool_metabolite_transport_rxn, transportReactionEquation); 

                                poolMetaboliteExchangeReactionEquation = strcat('EX_',strrep(listOfAnalyzedMetabolites{z,1},'[c]','(e)')); 
                                exchangeReactionEquation = sprintf('%s %s',extracellulatrPoolMetabolite, '<=>' );
                                currentModel = addReaction(currentModel, poolMetaboliteExchangeReactionEquation,exchangeReactionEquation); 

                                index_of_pool_extracellular_metaboilite = find(strcmp(extracellulatrPoolMetabolite, currentModel.mets),1);
                             
                                keggID = regexp(currentMetaboliteStructure{1}.ENTRY, 'C.....', 'match');
                                currentModel.metFormulas{index_of_pool_extracellular_metaboilite} = currentMetaboliteStructure{1}.FORMULA;
                                
                                if isempty(keggID)
                                    
                                else
                                    currentModel.metKEGGID{index_of_pool_extracellular_metaboilite} = keggID{1};
                                end
                                pool_met =  regexprep(regexprep(listOfAnalyzedMetabolites{z,1}, '[[c]]',''),'added_','');
                   
                                if ~strcmp(pool_met,elementsToExclude)
                                     listOfPoolMetabolites = [listOfPoolMetabolites; pool_met];
                                end  
                            end
                            
                            poolMetaboliteIndex = find(strcmp(listOfAnalyzedMetabolites{z,1}, currentModel.mets),1);
                            currentModel.metFormulas{poolMetaboliteIndex} = currentMetaboliteStructure{1}.FORMULA;
                            keggID = regexp(currentMetaboliteStructure{1}.ENTRY, 'C.....', 'match');
                            if ~isempty(keggID)
                                currentModel.metKEGGID{poolMetaboliteIndex} = keggID{1};
                            end
                        end
                    end
                end
            end

           if ~iFExist
               f = ['not exists ' pathway(y,:)]
                addedReactionIndex = find(strcmp(pathway(y,:), currentModel.rxns),1);
                [unbalancedRxns] = checkMassChargeBalance(currentModel,addedReactionIndex);

                if unbalancedRxns > 0
                    pathwayStatistics.UNBALANCED_REACT_COUNT = pathwayStatistics.UNBALANCED_REACT_COUNT + 1;
                end       
           end
        elseif regexp(pathway(y,:), 'C') == 1
            allMetabolitesInPathway = [allMetabolitesInPathway; pathway(y,:)];
        end
    end     
    
     if isForProduct
      %  objectiveRIndex = find(currentModel.c,1);
        addedRIndex = find(strcmp(currentModel.rxns, sourceTransportReaction));
        originalObjectiveReaction = evalin('base', 'originalObjectiveReaction');
sourceTransportReaction
       % currentModel = changeObjective(currentModel, sourceTransportReaction); %ne nado menjatj, budet originaljnoe izprediduwej modeli
        originalSolution = optimizeCbModel(currentModel); 

        pathwayStatistics.FBA_RESULT = originalSolution.f;
        pathwayStatistics.FBA_RESULT_OBJECTIVE_REACTION = originalSolution.x(addedRIndex);
        pathwayStatistics.FBA_BIOMASS_RESULT = originalSolution.x(originalObjectiveReaction);

        if  originalSolution.f == 0 || ...
            originalFlux >= originalSolution.f || ...
            originalSolution.x(originalObjectiveReaction) == 0 || ...
            originalSolution.x(addedRIndex) == 0
            zeroSolution = true;
        end

        
     else
        optimisationResult = optimizeCbModel(model);
        minimumBiomassFlux = optimisationResult.f * str2double(defines.minBiomassFluxProportion);

        objectiveRIndex = find(currentModel.c,1);
        addedRIndex = find(strcmp(currentModel.rxns,sourceTransportReaction));

        existingObjectiveReaction = currentModel.rxns(objectiveRIndex);
        model = changeRxnBounds(model, existingObjectiveReaction, minimumBiomassFlux,'l');
        currentModel = changeRxnBounds(currentModel, existingObjectiveReaction, minimumBiomassFlux,'l');

        if ~skipObjectiveReactionFluxAnalysis
            %original model - getting objective reaction flux
            model = changeObjective(model, model.rxns(objectiveReactionIndex - 1));
            currentModel = changeObjective(currentModel, currentModel.rxns(objectiveReactionIndex - 1));

            originalSolution = optimizeCbModel(model, optimizationObjective); 
            originalObjectiveValue = abs(originalSolution.f);

            calculatedSolution = optimizeCbModel(currentModel, optimizationObjective); 
            calculatedObjectiveValue = abs(calculatedSolution.f);

            pathwayStatistics.FBA_RESULT_OBJECTIVE_REACTION = calculatedObjectiveValue;
            pathwayStatistics.FBA_RESULT = calculatedSolution.x(addedRIndex);
            pathwayStatistics.FBA_BIOMASS_RESULT = calculatedSolution.x(objectiveRIndex);

            if  calculatedObjectiveValue == 0 || ...
                calculatedSolution.x(addedRIndex) == 0 || ...
                calculatedSolution.x(objectiveRIndex) == 0
                zeroSolution = true;
            end

            if calculatedObjectiveValue <= originalObjectiveValue
                zeroSolution = true;
            end
        else 
            currentModel = changeObjective(currentModel, currentModel.rxns(addedRIndex));
            calculatedSolution = optimizeCbModel(currentModel, optimizationObjective); 
            currentModel = changeRxnBounds(currentModel, sourceTransportReaction, calculatedSolution.f,'l');

            pathwayStatistics.FBA_RESULT_OBJECTIVE_REACTION = 0;
            pathwayStatistics.FBA_RESULT = calculatedSolution.f;
            pathwayStatistics.FBA_BIOMASS_RESULT = calculatedSolution.x(objectiveRIndex);

            if  pathwayStatistics.FBA_RESULT == 0 || ...
                pathwayStatistics.FBA_BIOMASS_RESULT == 0
                zeroSolution = true;
            end
        end
     end

    pathwayStatistics.MET_COUNT = pathwayStatistics.MET_COUNT - (((length(pathway(:,1)) - 1) / 2) - 1);
    pathwayStatistics.REACTIONS_WITH_KNOWN_THERMODINAMICS = reactionsNumberWithThermodinamics;
    
    
    %TODO CHECK FORMULAS AND IMPLEMENT COMBINED CALCLULATION
    reactionsWeight = calculate_reactions_weight( ...
        pathwayStatistics.REAC_COUNT, ...
        pathwayStatistics.EXISTING_REACT_COUNT, ...
        pathwayStatistics.UNBALANCED_REACT_COUNT, ...
        pathwayStatistics.REACTIONS_WITH_KNOWN_THERMODINAMICS);
   
    metabolitesWeight = calculate_metabolites_weight(...
        allMetabolitesInPathway, ...
        currentModel, ...
        pathwayStatistics,...
        listOfPoolMetabolites);

    listOfAllWeights = [pathwayStatistics.REAC_COUNT, ...
                        round(pathwayStatistics.FBA_RESULT*100)/100, ...
                        round(pathwayStatistics.FBA_BIOMASS_RESULT*100)/100, ...
                        round(pathwayStatistics.FBA_RESULT_OBJECTIVE_REACTION*100)/100, ... 
                        round(reactionsWeight * 100)/100, ...
                        metabolitesWeight];
    
                    
    if shouldModelBeReturned
        modelToReturn = currentModel;
    end
   % zeroSolution = false;
end



